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The thermodynamics of a charge-asymmetric lattice gas of positive ions carrying 
charge q and negative ions with charge —zq is investigated using Debye-Hiickel the- 
ory. Explicit analytic and numerical calculations, which take into account the for- 
mation of neutral and charged clusters and cluster solvation by the residual ions, 
are performed for z = 2, 3 and 4. As charge asymmetry increases, the predicted 
critical point shifts to lower temperatures and higher densities. This trend agrees 
well with the results from recent Monte Carlo simulations for continuum charge- 
asymmetric hard-sphere ionic fluids and with the corresponding predictions from 
continuum Debye-Hiickel theory. 
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I. INTRODUCTION 

The nature of critical phenomena in ionic systems has been a subject of numerous theo- 
retical studies in recent years.Bii'S Due to the long-range nature of Coulombic interactions, 
construction of a full renormalization group treatment, which was so successful in describ- 
ing critical behavior of non-ionic fluids, meets both conceptual and technical difficulties.0 
However, in recent years some progress has been achieved in obtaining physically reason- 
able, well-based mean-field theories for ionic systems-Si'llli These theoretical studies have 
been supported and, sometimes in substantial part, initiated by intensive Monte Carlo 
simulationji0i'illl0 of charged systems. 

To investigate the thermodynamics of ionic fluids, two main mean-field approaches have 
been developed. The first onellli extends the pioneering work of Debye and Hiickel0 (DH) on 
dilute solutions of strong electrolytes, while the second approach! relies on integral equations 
for correlation functions. Analysis of thermodynamic energy bounds and comparison with 
the best Monte Carlo estimates for the critical parameters suggests that the DH-based 
theory gives a better description of the thermodynamics of electrolytes, at least in the 
critical region.iHli 

The simplest model of ionic fluids, the so-called restrictive primitive model (RPM), con- 
siders a system of spherical equisized charged particles, half of them carrying a charge q 
and the other half with charge — g. The charge symmetry of this model plays a crucial role 
in the determination of its universality class and in the ability to obtain analytic solutions. 
This raises the question of how the breaking of the symmetry will affect the thermodynam- 
ics and critical properties of electrolyte systems. An important extension of the RPM is 
the charge-asymmetric primitive model, where the sizes of negative and positive particles 
are the same while absolute values of charges for positive and negative ions are different. 
Recent Monte Carlo simulationsi'0 have revealed that, as charge asymmetry increases, the 
critical temperature of the gas-liquid transition decreases, while the critical density Pc 
grows. However, most of the current theories give different predictions. Simple DH theory 
and the mean-spherical approximation (MSA) both predict no dependence on the asym- 
metry parameter In symmetric Poisson-Boltzmann and modified Poisson-Boltzmann 
integral equation methodsEl the charge asymmetry hardly changes Tc , while pc increases. 
However, the absolute values of the critical parameters are very different from Monte Carlo 
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estimates. A field theoretical approach by Netz and OrlandS predicts large increases in 
critical temperatures and, similarly, a large decrease in critical densities, in sharp contrast 
with computer simulations trends. The only theory that produces reasonable results for the 
effect of charge asymmetry on thermodynamics and criticality of ionic systems, as judged 
by comparison with Monte Carlo computer simulations, is the DH approach augmented by 
Bjerrum cluster formation and cluster- ion interactions (DHBjCI).0 

Lattice models, such as the Ising model, have played an important role in understanding 
criticality in non-ionic systems. In recent years, lattice models have also attracted the 
attention of researchers as a tool for investigating thermodynamics and criticality in Coulomb 
systems.0llllll0'0'il A systematic study of electrolytes on lattices, which utilizes the Debye- 
Hiickel approach, has been presented recently.0 In this work the thermodynamics of a d- 
dimensional system of equal numbers of positive and negative ions, i.e., a charge-symmetric 
lattice RPM, has been investigated. Specific calculations for Coulomb systems on three- 
dimensional lattices, which included ion pairing and ion-dipole interactions, predicted a gas- 
liquid phase separation at low densities. However, taking into account the lattice symmetry 
yielded a different scenario — the phase diagrams of electrolytes on simple cubic and body- 
centered cubic lattices show order- disorder phase transitions with a tricritical point, while 
gas-liquid phase separation is suppressed. The introduction of charge asymmetry in lattice 
models of ionic fluids tends to suppress the possibility of order-disorder phase transitions, and 
gas-liquid phase coexistence reappears, although the position of critical point may change. 

In this paper, we present a thermodynamic investigation of charge-asymmetric lattice 
models of electrolytes. By explicitly including the clustering of oppositely charged particles 
and ion-cluster interactions, we obtain phase diagrams for 2:1 and 3:1 lattice electrolytes, and 
we locate the critical point for the 4:1 ionic system. Our results accord well with the trend 
oMa.ed . .cen. Monte Ca* .n,uMio,. and the continnun, DHBjCI t,.eo..0 The pape. 
is organized as follows. In Sec. II we present an overview of our thermodynamic approach to 
multicomponent charged species mixtures and we outline the pure Debye-Hiickel theory. The 
full theory, which accounts for charged and neutral cluster formation and their interactions 
with the residual ions, for a 2 : 1 system is presented in detail in Sec. III. Sec. IV describes 
the general scheme of thermodynamic calculations for 3 : 1 and 4 : 1 lattice electrolytes. 
Finally, the results are discussed in Sec.V and our conclusions are given in Sec. VI. 
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II. LATTICE DEBYE-HUCKEL THEORY OF CHARGE- ASYMMETRIC 

ELECTROLYTES 

A. Thermodynamic Overview 

Consider a system of charged particles on a three-dimensional simple cubic lattice with 
a unit cell length a, which initially has A^^o ions carrying a charge —zq and zNq ions with 
a charge g, i.e., the total number of ions is = {z + l)iVo. Because of the electrostatic 
interactions, ions with opposite charges tend to form clusters. As a result, there will be many 
species present in the system: dimers, trimers, etc., with respective charges (— 2; + l)g, {—z + 
2)q, .., 0. If the number of particles of type i is given by Ni, then we define pi = Ni/V and 
Pi = Pi^o to be the number density and the reduced density of the i-th species, with vq = 
being the unit lattice cell volume. 

The Helmholtz free energy is central for the determination of the thermodynamic behavior 
of charge-asymmetric lattice electrolytes. It can be approximated by summing consecutively 
the free energies describing the interactions between different species, 

F = F''' + J2Fi, (1) 

i 

where F^"^ is the ideal lattice gas (entropic) term and Fi is the electrostatic energy of the 
i-th species. Once the reduced free energy density / = —F/kBTV is known, the reduced 
chemical potentials for every component fli = ii/ksT can be computed via 

/i, = -df/dp,. (2) 

Finally, the reduced pressure is given by 

p = p/kBT = f + ^PiPi. (3) 

i 

Then the possible phase equilibria are defined by matching pressures and chemical potentials 
for each component in different phases. 

In multicomponent systems with charged particles it is the electrochemical potential that 
must be equal in coexisting vapor (v) and liquid (/) phases,0 namely. 



Pi,v + qi(j)v = Pi,i + qi(f)h 



(4) 
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where is the electrostatic potential in the corresponding phase, where, in general, there 
is a nonzero Galvani potential difference A0 = (p^ — (pi between the phases.B However, for 
calculating phase equilibria in multicomponent systems of charged particles, it is more con- 
venient to use the single- component thermodynamic picture.0 Since every thermodynamic 
phase is electroneutral, the multicomponent system with N = [z + 1)Nq ions can be viewed 
as a single- component system of A^^o molecules, each of them consisting of one negative ion 
and z positive ions. Then phase equilibrium between the liquid and vapor at temperature 
T is ensured by 

P{T,p^) = P{T,pi) and ij{T , p^) = n{T , pi) , (5) 

where p„ and pi are the overall particle densities in gas and liquid phases, respectively, while 
p. = fi- + zp,+ . The pressure in each phase can still be calculated using Eq. @. This 
approach accounts for the electroneutrality of each phase, and utilizes only one chemical 
potential, which significantly simplifies calculations of phase diagrams. 



B. Pure DH theory of charge-asymmetric lattice electrolytes 

As a first approximation, assume that there is no clustering between oppositely charged 
particles, and that only free ions are present in the system. The free energy density can be 
written as / = f^'^ + f^^, where the first term is the entropic ideal gas contribution, which 
is given by 

f-id = _Pli^p* _ 1^ in(l - p*). (6) 

Vo Vo Vo 

The subscripts "+" and "-" denote positive and negative ions, respectively. Owing to overall 
electroneutrality in the system, the densities of free ions are related to each other by 

Pl = zp*/il + z), p*_=p*/{l + z). (7) 

The second term in the free energy density is the electrostatic contribution f^^ = 
/+(p+) + f-{p-), which is the result of ion-ion Coulombic interactions. By solving the lat- 
tice version of the linearized Poisson-Boltzmann equation for the electrostatic potential, and 
subsequently applying the Debye charging procedure, as was done for the charge-symmetric 
lattice model of electrolytes ,0 it can be easily shown that 



DH 
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d{x'), (8) 
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where 
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1 — I (cos ki + cos k2 + cos ^3) 



(9) 



is the integrated lattice Green's function for the simple cubic Iatticep3c3 and = k^o^ with 
the reciprocal squared Debye screening length given by 



Now the chemical potential for each type of ion can be easily calculated through Eq. 
which yields 



(10) 
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where we have defined the reduced temperature by 

DkBTa 



(11) 
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which is related to the reduced density as p* = x^T* /Att. Then for the chemical potential 
of the neutral "molecule," consisting of one negative and z positive ions, we obtain 
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(13) 



l-p*J ' " ' 3T* 

where the density-independent constant is C{z) = [zln z ~ {1 + z)ln{l + z)]. Having the 
chemical potentials and free energy, we can calculate the pressure using Eq. (^ to obtain 



pvo = -ln(l -p*) + — 



x'P 



6 



x"^ + Q 



P 



6 



+ 6 



rf(x2) 



(14) 



Note that the expression for the pressure is independent of the charge asymmetry, and the 
chemical potential of the "molecule" is {1 + z)/2 times the value for the 1:1 electrolytes 
(if we neglect the constant C{z) which does not affect the thermodynamics of the system). 
Then, the predicted phase separation and the critical point. 



T* = 0.1018, p* = 0.0996, 



(15) 



are the same as for the charge-symmetric lattice electrolytes. This result fully agrees with 
the pure DH theory in continuum space,0 which also predicts no change in the critical 
parameters for charge-asymmetric electrolytes. 
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III. LATTICE DEBYE-HUCKEL THEORY WITH CLUSTERING AND 
CLUSTER-IONS INTERACTIONS FOR 2:1 ELECTROLYTES 

A. Bjerrum clustering 

One of the main deficiencies of the pure DH theory, which describes the system of free 
positive and negative ions, is the total neglect of the ion clustering. Oppositely charged 
particles attract each other and form clusters in order to reduce the free energy. This process 
significantly decreases the number of free ions. In charge- symmetric lattice electrolytes in 
the critical region the number of neutral ion pairs is 2-3 times larger than the number of 
free ions,@ and one can expect that for charge-asymmetric ionic systems this ratio should 
be even larger because of the larger free energy gain for clusters with the increase of the 
parameter z. 

The importance of ion pairing in charged particles systems was first recognized by 
Bjerrum. @ In his original approach a cut-off distance between two oppositely charged ions 
was introduced to define a bound pair. Later, the definition of a pair became a subject 
of many discussions.B'i However, as was shown in a careful analysis by Levin and Fisher,!! 
the precise value of this cut-off distance has little influence on the critical parameters and 
coexistence curves (less than 0.5%). Meanwhile, the further problem of calculation of the 
electrostatic energy of the dipole particle poses another, new technical difficulty, since the 
ion pair does not possess spherical symmetry and the problem cannot be solved exactly 
analytically or numerically. By approximating the ideal bispherical exclusion zone by a 
symmetrically centered sphere. Levin and Fisheri succeeded in obtaining a precise, numeri- 
cally tractable solution. For three- or four- particle clusters in charge-asymmetric continuum 
ionic fluids, the same strategy also yields reasonable resultsjll although the technical com- 
plexity of computations increases signiflcantly. 

In lattice systems the situation is intrinsically simpler because of the discrete rather than 
continuous symmetry. In 1 : 1 electrolytes it allows one to deflne clearly an ion pair as 
two oppositely charged particles sitting on neighboring lattice sites.0 Similarly, for charge- 
asymmetric lattice electrolytes we can easily deflne different cluster conflgurations. These 
particles may be viewed as independent chemical species, and the processes of clustering can 
be considered as a set of chemical reactions between them. 
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To illustrate our approach, consider a 2 : 1 lattice ionic fluid in which clusters are allowed 
to form. Then, following Ref.0, in addition to the free positive particles with charge +g and 
free negative ions with charge — 2g, we suppose it suffices to consider three basic types of 
clusters: specifically we include dimers carrying a charge —q (with the number density pg) 
and two possible configurations of neutral trimers, linear and angular, with number densities 
Pg^ and pg^, respectively: see Fig.l. 

In the spirit of Bjerrum's approach, we assume first that neutral clusters do not interact 
with charged particles. Then these neutral particles contribute only to the ideal free energy, 
and the total free energy density is given by 

/ = f ^(p+, p_, P2, P3a, P3.) + /f + + /f . (16) 

To calculate the free energy of the charged dimers /^', we adopt a view of them as a 
combination of neutral symmetric (+, — ) dipoles, which occupy two neighboring sites, and 
single-site particles (with the charge —q) which sit on the top of negative part of the dipoles, 
i.e., (+,2—) = (+, — ) + (— ). Then the electrostatic free energy densities for the free ions 
and charged dimers are given by the corresponding expressions from the pure DH theory, 
although with a different inverse squared Debye screening length, namely, 

^' = ^{p^<f + ^P-<f+P2q% (17) 

To determine the ideal gas (entropic) contribution to the free energy one needs to know 
the corresponding densities and chemical potentials. As shown in Fig.l, we identify five 
different species in the 2 : 1 lattice Coulombic system. Let us define a as a set of free 
positive and negative ions, (3 as a. set of negatively charged dimers and free positive ions, 3a 
as the linear neutral trimers, and 36 as the neutral angular trimers. Then the clustering in 
these system can be described by three chemical reactions 

a ^ /?, (18) 
/3 3a, (19) 
/3 36, (20) 

with the corresponding equilibrium constants Ki,K2a and K2b- The chemical equilibrium in 
the system is described by the following relations between chemical potentials 

2p+ + p_ = p+ + P2 = P3a = Psb- (21) 
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The chemical potential for each particle can be computed by first separating the entropic 
and electrostatic parts, i.e., fii = /if'^ + /if, where the latter is determined simply by /if = 
—df^^/dpi- For the ideal-gas parts of the chemical potentials the situation is more complex. 
For free single-site ions the chemical potentials can be easily found applying the potential 
distribution theoremS or from simple entropic considerations, 

Ji'^ = Inp; - ln(l -p\- p*_ - 2p\ - 3p*, - Sp^J, (22) 
-p'^ = Inpl - ln(l -pl- p*_ - 2p; - Spl, - 3p;,). (23) 

However, for dimers and trimers there are no similar exact expressions. Nevertheless, these 
entropic contributions to the chemical potentials can be estimated by applying the Bethe 
approximation,!! which has been successful in the DH theory of 1 : 1 lattice electrolytes.0 
Note, that in lattice electrolyte systems phase transitions predominantly take place in the 
low-density regimes, and thus the error in using the Bethe approximation is expected to be 
small. 

Defining the activities of the particles via Zi, for every species we have 

z, = 0/Are'^S (24) 

where n is the number of sites occupied by the particle, Q is the corresponding internal 
partition function and A, is the de Broglie wavelength (see Ref.i). For the latter the equality 
holds A_,_ = A_ = A2 = A^a = Asfe. We refer to the Appendix A for a detailed calculation 
of the dimers and trimers activity using the Bethe approximation method, and we give here 
only the final expressions for activities, 

- f-'^J"!, „.„ e^- (25) 
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(l-p;-P--2p^-3pL-3pSj3' ' ^''^ 

= ^36(1 - \p*2f (,-yi.^-Pz^)(,f^^i (27) 

64(l-p;-p*_-2p^-3pL-3p*,)3 • V ; 

Once the activities are known, the chemical potentials are derived by utilizing Eq.(^) for 
each species, and the free energy is obtained by integrating the chemical potentials. Finally, 
we arrive at the chemical potential of an a- "molecule" (which consists of two positive and 
one negative ions) in the form 

pf = Inp; + 2 lnp*_ - 3 ln(l - p^ - p*_ - 2pl - 3p*, - 3p*,) + (28) 
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and the reduced pressure 



pD^^H, = - ln(l -pl-p*_- 2p; - 3p*„ - 3p* J + 3 ln(l - p;/3) + 3 ln(l - plJ2) 



+ 31n(l-py2) 



(29) 



2-P3a 

The electrostatic part of the chemical potential and the pressure have the same form as in 
the pure DH theory (p!4D, namely, 
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(30) 
(31) 



where x^ = K^a^ with given by (0). 

In terms of activities, the chemical equilibrium conditions (^) can be presented in the 
form of laws of mass action 



Z2Z+K2b 



Z2Z+, 
Z3b, 



(32) 
(33) 
(34) 



As can be seen from the chemical equilibrium conditions (|2l|) and the definition of activities 



(1^), the association constants in (|3^)- (|3^) are equal to the internal partition functions of 
dimers and trimers. 



The above definitions of the association constants lead to the expressions 

n 



(35) 



(36) 



where Uji = qjipj is the electrostatic energy of jth ion of a cluster particle i in the potential 
yjj, which is due to the other ions entering the multimer cluster. The coefficient Cj in Eq.(p6D 
is an entropic factor which takes into account all allowable configurations of the cluster. The 
electrostatic potentials (pj can be determined by solving the lattice version of the linearized 
Poisson-Boltzmann equation 



(37) 
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in which and are the charge and the position of the kth. ion entering the multimer 
cluster. Since lattice Coulomb potentials can be calculated exactly in numerical terms,^ we 
obtain the following expressions for association constants 

^1.08152" 



K2b 



6f exp 
vq exp 
4t>o exp 



0.812033 

p* 

"0.734737 



(38) 
(39) 
(40) 



Substituting 
chemical equilibrium 



and the expressions for activities (P5|)-(P7D into the expressions for 
2D-f0) yields a set of equations, which define implicitly the dimer 
and trimer densities Psa ^"^^ Psb terms of the monomer densities and p*_. The 
electroneutrality of the system requires that p^ = 2p!_ + p2. Then the chemical potential 
and the pressure can be expressed in terms of only one variable, the total reduced density 
P* = P+ + P- + 2p2 + 3p3a + 3p3^, which allows for the construction of the coexistence curve: 
see Fig. 2. 

As in the case of charge-symmetric continuumi and latticeH electrolytes, the coexistence 
curve in the DHBj approximation has an unphysical banana-like shape. This is related to 
the fact that, as the temperature is lowered, the number of neutral clusters quickly grows, 
and this depletes the number of free charges present in the system. Since it is the density of 
free charges that plays the role of the order parameter and governs the gas-liquid transition, 
the phase separation takes place at higher overall densities. At this level of approximation, 
neutral clusters are electrically inactive, and hence contribute only to the hard-core part of 
the free energy, which is the same for both phases. Thus, DHBj theory simply superimposes 
the pressure of ideal gas of clusters on the pure DH pressure, and both sides of the coexistence 
curve shift to higher densities by equal amounts. The critical density is now substantially 
higher, p* ^ 0.0807, while the critical temperature slightly decreases, T* ^ 0.099. 



B. Cluster-Ion interactions 



The predictions of DHBj theory for lattice ionic systems are thermodynamically unreason- 
able and should be corrected by taking into account the effects of interactions between multi- 
mers in clusters and free ions.iil As shown earlier for continuum and lattice electrolytes,liii 
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these solvation effects eliminate the unphysical banana-shape phase coexistence curves. It is 
reasonable to expect that these cluster-ion interactions are even more important in charge- 
asymmetric ionic systems due to a larger fraction of neutral clusters in equilibrium with free 
ions. 

The exact calculations of interactions between the multimers and free ions are very com- 
plicated. Instead we use a reasonable assumption to obtain closed analytic expressions. 
We approximate the neutral clusters as a set of overlapping non-interacting symmetric 
dipoles. For example, the neutral trimers can be viewed as a combination of two dipoles, 
i.e., (+, 2— ,+) = (+, — ) + ( — ,+). The charged clusters, as we discussed above, are ap- 
proximated as symmetric dipoles overlapping with free ions, i.e., (+,2—) = (+, — ) + (— ). 
We expect that differences between our approximation and the exact cluster-ion interaction 
contribution to the free energy to be small since the corrections correspond to higher order 
dipole-dipole interactions, and thus can be neglected. 

According to our approximation the reduced density of dipoles from all cluster particles 
in the system is equal P2 + 2p3^ + 2p*^^^. Since the energy of solvation of a single free dipole 
in the lattice electrolyte system is known exactly,^ the cluster-ion interactions contribute 
to the free energy density as 



/^^ = (p; + 2p*„ + 2p: 



where 



^^'2lDkBTvl 
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+ 6 
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(41) 
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Then the corresponding contributions to the chemical potentials and pressure are given by 



27r2(p* + 2p*, + 2py 
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(43) 



(44) 



/i^^x^T747r, 

which must be added to the values provided by Eqns. (p8|) and (pQ]) in order to obtain the 
complete expressions for the chemical potential and pressure. 

The full DHBjCI theory predicts a phase separation as exhibited in Fig. 2. As in the 
charge-symmetric lattice and continuum electrolytes,iS taking into account the solvation 
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of clusters by the residual free ions eliminates the unphysical banana-like shape of the co- 
existence curve. The critical parameters are now given by 

T* = 0.08735, p* = 0.05001, = Pc/ PcksT = 0.2431. (45) 

Comparison with the critical parameters of 1:1 lattice ionic system indicates that the critical 
temperature is about 11% lower, while the critical density is 1.67 times higher. These trends, 
the decrease in the critical temperature and inrease in the critical density, are in agreement 
with the results of continuum calculations for charge-asymmetric electrolytes.lll 

It is interesting to estimate the relative amounts of different species in the critical region. 
In terms of molar fractions yi = nip*/p* [rii is the size of the particle), our theoretical 
approach predicts 

y+ = 0.1002, y_ = 0.0388, y2 = 0.0452, (46) 
ysa = 0.0075, and y^b = 0.8083. 

As expected, the fraction of neutral trimers significantly exceeds the fraction of free charges. 
Also, the number of charged dimers is very low, which can be attributed to their propensity 
to combine with free ions to form energetically more favorable neutral trimers. What is 
surprising, at first glance, is that the ratio of linear trimers to angular trimers is much 
less than unity. Indeed, because electrostatic interactions favor a linear arrangement of the 
ions in a trimer cluster, these clusters are more energetically stable and should prevail over 
angular trimers. However, these naive arguments do not take into account the entropic 
considerations. First of all, the angular trimers have more different possible configurations 
than the linear trimers. In addition, they are more compact and can be packed more densely. 
As a consequence, there are more ways to arrange angular trimers on the lattice and they 
dominate trimer clusters. 

IV. LATTICE DHBJCI THEORY FOR 3:1 AND 4:1 ELECTROLYTES 

The thermodynamic calculations for z = 3 asymmetric lattice electrolytes can be per- 
formed following the method presented in full detail in Sec. III. As shown in Fig. 3, there are 
seven basic clusters to be considered in this system, namely, single-site positive and negative 
ions, dimers with the charge —2q, linear and angular trimers (with the charge — g), and two 
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types of neutral tetramers, which are connected by the network of six chemical reactions. 
Note that we do not consider clusters with bonds between the same charges. For example, we 
assume that charged trimers have the configurations (+,3— ,+), but not (+,+,3—). Indeed, 
these particles are less stable and their contributions to free energies can be neglected. 

In addition, as in the case of 2:1 lattice Coulombic systems, we view the clusters as com- 
binations of noninteracting (+,— ) symmetric dipoles and single-site charges. The activities 
of multimers are calculated by applying the Bethe approximation method outlined in the 
Appendix. These approximations allow us to calculate the free energy of the system as a 
sum of the corresponding entropic, electrostatic and cluster solvation terms for each particle. 
The resulting phase coexistence curve is given in Fig. 4. The critical parameters are 

T* = 0.0688, pI = 0.0847. (47) 

The molar fractions of different clusters in the critical region are 

y+ = 0.08872, y_ = 0.00347, y2 = 0.0088, (48) 
l/3a = 0.0033, 1/36 = 0.2049, y^a = 0.3043, y^^ = 0.3864. 

As expected, neutral clusters again dominate. 

Analogous thermodynamic calculations can be done for 4:1 lattice electrolj^es. However, 
the number of different types of clusters and the number of chemical reactions between them 
become fairly large, and the full thermodynamic analysis is difficult to complete. Instead, 
we focus on the critical region of the system where calculations can be completed. The 
resulting critical parameters are 

T* = 0.060, pI = 0.148. (49) 

The relative density of all neutral clusters here is about 87%. 

Full phase diagrams for 2; : 1 lattice ionic systems (^=1,2 or 3) calculated using DHBjCI 
approach are exhibited in Fig.4. 

V. DISCUSSION 

Our analysis of charge-asymmetric lattice electrolytes using Debye-Hiickel theory with 
Bjerrum clustering and cluster-ion interactions indicates that charge asymmetry strongly 
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influences the thermodynamics, especially in the critical region. It is found that critical 
temperatures decrease, while critical densities increase with charge asymmetry. Our theo- 
retical predictions for critical parameters are compared in Fig. 5 with the results of continuum 
Monte Carlo simulationsi00 and with the predictions of the continuum DHBjCI theory.llZl 

Clearly, the overall consistency between the lattice and continuum DHBjCI approaches 
and computer simulations results indicates that the Debye-Hiickel method correctly captures 
the physics of phase transitions in these ionic systems. In Coulomb systems with charge 
asymmetry the formation of clusters is strongly favored. Clusters contribute to electrostatic 
interactions much less than free ions, and the effective electrostatic energy per particle 
decreases. However, the phase separation is driven by charged particles, and the temperature 
is normalized by the strongest electrostatic interaction between the positive and negative 
ions: see Eg. ([T^) . Therefore, the reduced critical temperature falls. 

Cluster formation also significantly depletes the number of free ions. At the critical point, 
the molar fractions of free single-site positive and negative ions are 21.0%, 13.8%, 9.2% and 
4.2% for 1:1, 2:1, 3:1 and 4:1 lattice electrolyte systems, respectively. As a consequence, 
larger overall densities are required to achieve the phase separation, and the reduced critical 
density increases. 

At the same time, the relative numbers of neutral particles as a function of charge asym- 
metry shows a non-monotonic behavior. The molar fractions of neutral clusters are equal to 
79% in 1:1 ionic system, 81% in 2:1 electrolytes, 69.1% in 3:1 ionic system and about 88% in 
4:1 electrolytes. The decrease of the molar fraction of neutral clusters in the critical region in 
3:1 system as compared to 2:1 electrolytes has also been found in continuum calculations.^ 
The actual number of neutral clusters depends on two factors: the free energy gain of cluster 
formation and the temperature (through the equilibrium constant). It is possible that in 3:1 
electrolytes the lowering of the critical temperature is not enough to overcome the repulsion 
between the positive ions in neutral tetramers, while in 4:1 ionic systems the decrease in 
critical temperature is enough to favor more strongly the formation of neutral clusters. 

The absolute values of the critical temperature in our charge-asymmetric lattice models is 
larger than the Monte Carlo simulation datai'0 by a factor of 1.68 — 1.94, with lesser factors 
corresponding to higher z. This agrees with the general feature of lattice models to over- 
estimate the values of critical temperatures in comparison with continuum models. Note, 
however, that the simulations are performed with the continuum potential, while our theory 
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employs the lattice Coulombic potentials. In addition, the reduced temperature (|I^) is 
normalized by the energy of two ions interaction via continuum or 1/r Coulombic potential. 
However, the lattice potential differs from 1/r at short distances,^ and normalizing the 
temperature by the energy of lattice Coulombic interactions of nearest neighbors would 
lower the temperatures by an additional 8%. 

The absolute values of the critical densities in our theoretical approach are only 38-69% 
of the current Monte Carlo estimates.i'0 This fact reflects the discrete nature of the lattices 
and that not all possible clusters have been taken into account. 



VI. CONCLUSIONS 



We have investigated the effects of charge asymmetry on thermodynamics and critical 
properties of the lattice ionic systems using Debye-Hiickel theory with Bjerrum clustering 
and cluster-ion interactions. Phase diagrams for z : 1 lattice electrolytes have been obtained 
for z = 2 and 3, while for z = 4 the location of the critical point was determined. Our results 
agree well with the Monte Carlo simulations and with the continuum Debye-Hiickel theory,llll 
predicting that the increase in the charge asymmetry lowers the critical temperature and 
increases the critical density. 

Our theoretical approach may be extended in several directions. In this paper we investi- 
gated charge asymmetric electrolyte systems on the simple cubic lattices. It is interesting to 
consider the lattices with different symmetry such as body-centered cubic and face-centered 
cubic lattices, for which thermodynamic calculations for charge symmetric systems are al- 
ready performed.^ Another interesting question is how the charge asymmetry will affect the 
ionic systems in two dimensions, where interesting phase transitions may appear. 
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Appendix A: THE BETHE APPROXIMATION FOR MULTIMERS 



Consider a system of multimers (each particle occupies m lattice sites), which interact 
with each other only via on-site exclusion, on a simple cubic lattice with the total number 
of sites A^. We can calculate the entropy of the multimer distribution following a direct 
method of counting probabilities in the Bethe approximation employed by Nagleii. (For 
other methods of calculating the partition function for lattice systems see Refs.0'i.) The 
Bethe approximation does not take into account any correlations around cycles and thus it 
is exact for lattices with only tree-like paths (Bethe lattices). However for SL* cubic lattices 
this approximation gives fairly good results. 

Consider, first, a system consisting only of dimers. Then the number of ways to arrange 
N2 = N p*2 dimers on the simple cubic lattice can be approximated byil 



W2{pl 



N 
2Npl 



2p| 
6 



1/2 



2p| 
6 



5/2' 



2Np% 



N(l~2pl) 



(Al) 



Here the first square bracket represents the total number of ways to arrange N p2 dimer 
vertices on the lattice, taking into account six possible orientations for each dimer, provided 
one its vertex is fixed. The second bracket gives the probability that a dimer vertex config- 
uration is compatible with all its nearest-neighbors' vertex configurations. Here 1p\l^ is the 
probability of finding dimer's edge along certain lattice direction starting from a given point, 
and (1 — 2p2/6)^ ensures there is only one dimer at this point. The third square bracket 
factor is the probability that no dimers occupy A^(l — 1p*^ empty lattice sites. The square 
roots are taken in order not to compute the probabilities twice. The Bethe approximation 
for dimer's activity on the sc lattice then yields 



exp 



1 d 



ini^2(p;; 



p;/3(l-p;/3) 
(l-2p*^2 



(A2) 



N dpi J K--^H2) 

Similar calculations can be performed for the systems containing trimers, tetramers or 
higher order multimers. However, since now not all multimer vertices are equivalent, we 
have to account for each type of vertex separately. 

In particular, consider the system with A3 = Np^ trimers. Then there are 3Np^ trimer 
vertices, with Np^ "center" and 2A^p3 "end" vertices. For linear trimers the number of 
possible orientations is 3 if we fix a "center" vertex and 6 if we fix an "end" vertex. Then 
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the number of ways to arrange linear trimers on the lattice is given by 



N 
3Npl 



2Np* 



X 



^6 J \ 6 



4/2' 



(A3) 



Again the first square bracket gives the maximum possible number of distinct trimer vertices 
configurations, the second square bracket accounts for the probabihty that all "end" vertex 
configurations are compatible with their nearest-neighbors' configurations, the third bracket 

corresponds to the "center" vertices' compatibility, and the last factor ensures that there 
are A'^(l — 8^3) lattice sites free of trimers. Then the activity of linear trimers is given by 

P3a(l - 3pL/6)\ -l/(2-.^„) 



(A4) 



For angular trimers, the only difference is that now there are 24 distinct trimer orientations 
if we start from the "end" vertices, and 12 orientations if we start from the "center" vertex. 
This yields 



Zsb — 



1 P3(.(i - 3py6)2 



-l/(2-p^6) 



64 (l-3p*J3 

Similarly, the activity of tetramers is given by the following expression 

, ^4(1 - '2pi/3y i/^i_2p*/3) 



(A5) 



(A6) 



where the coefficient C equals to 1/96 for planar tetramers and C — 1/128 for non-planar 
tetramers. 
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Figure captions. 

Fig.l. Different charged and neutral particles and chemical reactions in 2:1 lattice 
electrolyte systems. 

Fig.2. Phase diagrams for 2:1 lattice electrolyte: (a) pure DH theory, (b) DHBj 
approximation, (c) DHBj CI theory. 

Fig.3. Different charged and neutral particles and chemical reactions in 3:1 lattice 
electrolyte. 

Fig.4. Phase diagrams of (a) 1:1, (b) 2:1, (c) 3:1 lattice electrolytes in DHBjCI theory. 

Fig.5. A) Critical temperatures T* and B) critical densities p* as functions of charge 
asymmetry: (a) our lattice DHBjCI predictions, (b) continuum DHBjCI predictions,^ (c) 
Monte Carlo simulations 
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Fig.2 



24 




25 




Fig.4 
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Fig. 5 A 
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Fig.5B 



